Eddie module stream discharge from handout

This script examines the discharge in Neversink Cr.

This script uses the excel sheet provided in the Eddie module. It is possible to download data from USGS and put it in a new excel sheet following the same format and rename the parts of the script. You can also use the USGS retrieval package to download up to date data directly from USGS.

Install libraries.

These are the libraries that are used in the script. These are installed only once and then can be commented out so they do not run in future execuations of the script.

# https://github.com/USGS-R/dataRetrieval

# general R packages to install if you dont have them
# remove the "#" if you need to run the installations

# install.packages("devtools") # essential in installing other thigs
# install.packages("tidyverse") # installs a lot of things and ggplot
# install.packages("scales") # allows great scale formatting on ggplot
# install.packages("lubridate") # makes working with dates easier
# install.packages("janitor") # clean names of columns and other things
# install.packages("readxl") # read in excel files
# install.packages("plotly") # interactive plots
# install.packages("skimr")
# install.packages("patchwork") # multiple plots

# specific to this module
# install.packages("dataRetrieval") # USGS Data Retreiveal Method

Load Libraries

This code loads the libraries for working with the data

# load these every time  
library(tidyverse) # the tidyverse library with many sub libraries

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.2.1     ✓ purrr   0.3.3
✓ tibble  2.1.3     ✓ dplyr   0.8.3
✓ tidyr   1.0.0     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.4.0
── Conflicts ──────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(scales) # working with ggplot scales 

Attaching package: ‘scales’

The following object is masked from ‘package:purrr’:

    discard

The following object is masked from ‘package:readr’:

    col_factor

library(lubridate) # working with dates

Attaching package: ‘lubridate’

The following object is masked from ‘package:base’:

    date

library(janitor) # cleans up excel imports

Attaching package: ‘janitor’

The following objects are masked from ‘package:stats’:

    chisq.test, fisher.test

library(readxl) # read in excel files
library(skimr) # does summary stats and other things
library(plotly) # makes interactive plots

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(patchwork) # combines many plots together

# # load the USGS package
library(dataRetrieval) # utilities for downloading and working with USGS data

Source of Neversink data

This is the data that is in the excel sheet that is used in this script.
https://waterdata.usgs.gov/ny/nwis/uv/?site_no=01435000&PARAmeter_cd=00065,00060,63160

This is all of the data that is available for the site.
https://waterdata.usgs.gov/nwis/inventory/?site_no=01435000&agency_cd=USGS

Read in USGS Data

This allows you to read in the daily mean data for the USGS site.
Note Enter the following information and you can change this to get different sites or time ranges:
site_no: this is the site that you want to import
par_code: these are the paramter code numbers - 00060 is discharge
You can retreive other variables using:
c(“00060”, “00065”) 00065 gage height start.date:
end.date:

# The site is Neversink and the site number is 01435000
# USGS 01435000 NEVERSINK RIVER NEAR CLARYVILLE NY

site_no <-          "01435000"
par_code <-         "00060"
start.date <-      "1800-01-01"
end.date <-        "2050-01-01"

neversink_daily.df <- readNWISdv(siteNumbers = site_no,
                        parameterCd = par_code,
                        startDate = start.date,
                        endDate = end.date)
  
# rename the columns
neversink_daily.df <- renameNWISColumns(neversink_daily.df)

we need to make a date column to plot with

neversink_daily.df <- neversink_daily.df %>%
  select(agency_cd, 
         site_no, 
         date = Date, 
         discharge_cfs = Flow) 

we need to make a date column to plot with

neversink_daily.df <- neversink_daily.df %>%
  mutate(month = month(date))

Plot all data

q.plot <- neversink_daily.df %>% 
  ggplot(aes(date, discharge_cfs)) +
  geom_point(size=0.1) +
  geom_line()
q.plot

Plot February data

feb.plot <- neversink_daily.df %>% 
  filter(month == 2) %>%
  ggplot(aes(date, discharge_cfs)) +
  geom_point(size = 0.1) +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
feb.plot

NA

We can make it interactive

ggplotly(feb.plot)

August data

aug.plot <- neversink_daily.df %>%
  filter(month == 8) %>%
  ggplot(aes(date, discharge_cfs)) +
  geom_point() +
  geom_line() +
  labs(y="Mean Monthly Discharge m^3/sec", x="Date")
aug.plot

Compare the plots

feb_aug.plot <- feb.plot +
                aug.plot +
                plot_layout(ncol =1)
feb_aug.plot

add trend lines as straight lines

feb_aug.plot <- feb.plot + geom_smooth(method="lm") +
                aug.plot + geom_smooth(method="lm") +
                plot_layout(ncol =1)
feb_aug.plot

How to see regression statistics

feb.model <- neversink_daily.df %>% 
  filter(month == 2) %>%
  lm(discharge_cfs ~ date, data=.)
summary(feb.model)

Call:
lm(formula = discharge_cfs ~ date, data = .)

Residuals:
   Min     1Q Median     3Q    Max 
-129.4  -86.5  -55.4    0.3 5925.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.585e+02  5.282e+00  30.007   <2e-16 ***
date        1.393e-03  5.713e-04   2.439   0.0148 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 233.5 on 2258 degrees of freedom
Multiple R-squared:  0.002627,  Adjusted R-squared:  0.002186 
F-statistic: 5.948 on 1 and 2258 DF,  p-value: 0.01481

aug.model <- neversink_daily.df %>% 
  filter(month == 8) %>%
  lm(discharge_cfs ~ date, data=.)
summary(aug.model)

Call:
lm(formula = discharge_cfs ~ date, data = .)

Residuals:
   Min     1Q Median     3Q    Max 
 -91.5  -61.5  -36.9   -2.6 7104.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 8.495e+01  4.794e+00  17.720  < 2e-16 ***
date        1.987e-03  5.152e-04   3.856 0.000118 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 220.2 on 2477 degrees of freedom
Multiple R-squared:  0.005966,  Adjusted R-squared:  0.005564 
F-statistic: 14.87 on 1 and 2477 DF,  p-value: 0.0001184